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Abstract. - To synthesize proteins in a cell, an mRNA has to work with a finite pool of ribosomes. 
When this constraint is included in the modeling by a totally asymmetric simple exclusion process 
(TASEP), non-trivial consequences emerge. Here, we consider its effects on the power spectrum 
of the total occupancy, through Monte Carlo simulations and analytical methods. New features, 
such as dramatic suppressions at low frequencies, are discovered. We formulate a theory based 
on a linearized Langevin equation with discrete space and time. The good agreement between 
its predictions and simulation results provides some insight into the effects of finite resoures on a 
TASEP. 



Introduction. — Unlike systems in thermal equilib- 
rium, those driven far from equilibrium are poorly un- 
derstood. Even when they settle into time independent, 
steady states, there is no comprehensive framework that 
provides, say, the stationary distributions. Yet, such sys- 
tems are ubiquitous in nature, from small biological sys- 
tems like the cell to global networks like the internet. In 
order to gain some insight into the physics of such systems, 
we may first turn to simple models, from which more com- 
plex theories may be built. The totally asymmetric simple 
exclusion process (TASEP) [1] is a prime example, play- 
ing a role similar to the Ising model for the understanding 
of phase transitions in equilibrium statistical mechanics. 
Since its stationary distributions are known analytically, 
many interesting macroscopic properties can be computed 
exactly [2] . Beyond the simplest TASEP, many generaliza- 
tions have been proposed for modeling a variety of systems 
in science and engineering, from protein synthesis [3] and 
surface growth [4] to vehicular traffic [5] . 

Despite these advances, many of its behaviors continue 
to present us with surprises. One recent example is the 
power spectrum associated with the total particle occu- 
pancy showing dramatic oscillations (in the frequency do- 
main) [8]. Although the predictions from a linearized 
Langevin equation are found to agree well with simula- 
tion data, puzzles concerning some of the fit parameters 
linger [8]. In this study, we pursue this quantity further, 
but in a different setting. In modeling of protein synthesis, 
a TASEP represents a mRNA while the particles represent 



ribosomes that bind to one end of the mRNA, move to the 
other end, and detach. In solvable TASEPs, particles en- 
ter the system with a fixed rate, as if coupled to an infinite 
reservoir. Yet, in a real cell, the number of ribosomes is 
clearly finite. Thus, the effects of feedback from "ribo- 
some recycling" [10] deserve attention, especially if the 
effective binding rate may depend on the concentration of 
the available ribosomes. Recently, two studies investigated 
the effects of a finite pool of particles on TASEPs [11, 12], 
given reasonable assumptions on how the effective entry 
rate depends on the numbers in the pool. Here, we ad- 
dress another issue: How are the power spectra affected? 
Using Monte Carlo simulations, we found serious suppres- 
sions at low frequencies. In the remainder of this letter, 
we briefly discuss the model, the simulation results, and a 
theoretical explanation for this suppression. We end with 
an outlook for further research. 

Model. — The original open TASEP, which will be re- 
ferred to here as "ordinary," consists of a one-dimensional 
lattice with sites labeled by a; G [1, L], Each site may be 
occupied by at most one particle. In each time step, a ran- 
domly chosen particle hops to the next site (x — > x + 1), if 
the latter is empty. If chosen, the last particle leaves the 
lattice with probability /3, while the first site, if empty, 
will be filled by a particle with probability a. Thus, the 
total occupancy is a fluctuating quantity in time: N (t). 
In the steady state, the time average (N) and overall den- 
sity (p = N/L) settle to constants. In the thermody- 
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namic limit, three phases exist: (i) maximal current (MC) 
with p = 1/2 for a,/3 > 1/2, (ii) high density (HD) with 
p = 1—0 for j3 < 1/2 and a > f3, and (iii) low density (LD) 
with p — a for a < 1/2 and a < (3. At the HD-LD phase 
boundary (a = /3 < 1/2), a LD region (for smaller x) co- 
exists with a HD region, separated by an interface with 
microscopic width, known as the "shock" or the domain 
wall (DW) . As the DW performs a random walk through- 
out the lattice, the overall density suffers macroscopic fluc- 
tuations, but p is restored to 1/2 over long times. Though 
this HD-LD boundary is not a phase, it is often referred 
to as the "shock phase": (SP). 

In the constrained TASEP, we connect the lattice to 
a "pool" of N p particles. When the last particle leaves 
the lattice (with rate /3), it is recycled back into this pool, 
while a particle enters the lattice (if the first site is empty) 
with a A p -dependent rate a e //- Thus, the total number 
of particles in the system, N to t — N + N p , is a constant. 
Following previous studies [11, 12], we use 

(*eff (N p ) = atanh (N p /N*) (1) 

where N* marks a cross-over scale, chosen typically as 
pL. This form is motivated by the following. If there 
are few ribosomes available, the binding rate should be 
proportional to their concentration. Thus, ot e ff cx N p for 
small N p . On the other hand, regardless of how abundant 
the ribosomes are, the binding rate should not exceed some 
intrinsic value, which we modeled by a. For convenience, 
we will label the pool as a site (x = 0) on a periodic lattice 
with L + 1 sites. Clearly, there is no exclusion for this site 
(unrestricted N p ) and the rules of hopping into/out of it 
differ from the rest. 

For our simulation study, time t is measured in Monte 
Carlo Steps (MCS). In one MCS, L+l nearest-neighbor 
pairs of sites are chosen randomly for an update attempt. 
Starting with all particles being in the pool, the system is 
allowed to reach steady state by waiting 100k MCS typi- 
cally. Thereafter, the total number of particles on the lat- 
tice, N (t), is measured every t MCS until T MCS. Using 
mostly £,T = 10 2 , 10 6 , we compute the Fourier transform 
N(m) = Er=i N ( iT ) e 2mmTl / T with m E [1, T/l] . Typi- 
cally 100 such runs are carried out and our power spectrum 
is defined as I(m) = (|A(m)| 2 ), where (...) indicates the 
average over these runs. 

Simulation results. — With appropriate a and (3, the 
ordinary TASEP settles into a steady state: MC, LD, HD, 
or SP. For the constrained TASEP, we have an additional 
control parameter: N tot . In [11], the effects of tuning N tot 
on p (and the overall current) were reported. Not surpris- 
ingly, the TASEP is essentially in a LD phase when Ntot 
is small, regardless of (a,fl). As N to t is increased, more 
interesting crossover behaviors were found [11]. Here, we 
turn to /(to), which provides information on time correla- 
tions (N (t) N (t')). To distinguish the two power spectra, 
we will use I a for the ordinary TASEP [8] and I c for the 
constrained case here. In all cases, I c is severely smaller 
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Fig. 1: (Color online) Upper: Comparison of the power spectra 
of the ordinary and constrained TASEP with p = 0.1 and L = 
1000. Lower: Constrained power spectra for p — 0.1 with 
L = 1000 with different a e tt. The a' e ff = line corresponds 
to the ordinary TASEP with a constant a. 

than I Q at low frequencies. We will focus on two regimes: 
(i) the LD and (ii) crossing over to the HD, through the 
SP. 

For the LD regime, I {m) displays dramatic oscilla- 
tions which depend on L and a. Similar oscillations 
are found in I c (m), except for the suppression at low 
m's. An example is shown in fig. [T] where [3 = 0.9 and 
L = 1000 in both cases, while a = 0.1 for I (m) and 
(a,N tot ,N*) = (0.3,205,300) for I c (m). For large m's, 
the two spectra are approximately equal, since we have 
chosen the a for I Q to be approximately the average a e f / 
for I c . This behavior can be intuitively understood: On 
short time scales, fluctuations in N (t) do not have time 
to traverse the entire lattice and to contribute to the feed- 
back, via a e ff (Not — N). By contrast, fluctuations over 
longer times will affect ct e ff and so, 7 c (m) at low m. Fur- 
ther, we can expect a significant role from the "stiffness" 
of a,,ff, i.e., a' eff = d Np a e ff\^ , where N p = N tot - pL 
is the average pool occupation. To explore further how 
I c is affected by a range of cJ^p we investigated eq. |T]) 
with different N*'s (and other parameters so as to keep 
ct e ff close to 0.1). The results are illustrated in fig. [1] 
and confirm our belief that a larger ct' e ff induces a larger 
suppression. 

For the other regime - crossover through SP to HD, I c 
is also severely suppressed at low m, though the detailed 
structures are less interesting. We should emphasize that, 
in this SP-like regime, the average a e ff is equal to f), for a 
range of N to t values. Thus, we may expect I c ~ m~ 2 , just 
like in I Q [8]. Similar to above, this expectation is indeed 
valid, but only for large m. In the suppressed region at 
low m, I c approaches a constant value. 

Theoretical considerations. Since / (to) carries 
information on time correlations of N (t), even I a cannot 
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be accessed from the exact stationary distribution [2] . In- 
deed, finding I Q would be still a serious challenge [7], even 
if all the eigenvalues and the (left and right) eigenvectors of 
the Liouvillian (of the master equation governing TASEP) 
were explicitly known [6]! Thus, we will attempt to under- 
stand the suppression through a semi-phenomenological 
approach, similar to the one in [8]. Also based on a 
Langevin equation for the local particle density p(x,t), 
our approach here is improved over [8] in several ways. 
One is the use of discrete space and time points [9], thus 
avoiding any issues associated with UV cutoffs. Another 
concerns the use of a ring of L + 1 sites rather than open 
boundaries, so that any issues associated with boundary 
currents and noise are avoided. Finally, due to the con- 
straint, N (t) — N to t — N p (t), so that I (to) is directly 
related to the fluctuations of p(0, t), which is the pool. 

For simplicity, we define A t p(x, t) = p(x, t + 1) — p(x, t) 
and suppress t when no ambiguity exists. For all but the 
"exceptional sites," i.e., x = 0, 1, and L, we may write 

A t p(x) = p(x-l)[l-p(x)]-p(x)[l-p(x + l)} 

+£(SB-M) (2) 

where £(x, t) is the noise associated with the jump x — ► 
x + 1. Next, we define £ = £(L) — £(0) and write 

A t p(0) = (3p(L) - a eff (p(0)) [1 - p(l)} + 1 (3) 

and similar equations for p{\) and p(L). The conservation 
constraint, ^2 X A t p(x) = 0, can be verified. In principle, 
it is possible for p(x ^ 0) to exceed unity (due to the 
Gaussian noise and discrete time), but we will disregard 
this complication here. 

To simplify the computation, we choose a, /?, and N to t 
such that the system is in a LD state with a flat (average) 
profile, p. For fluctuations in the lattice (x € [1, £]), we 
write (p(x, t) — p(x, t)—p and, for the pool, £(i) = N p (t) — 
N p (due to the conserved dynamics, C + Sz f = *-*)• Next, 
we expand the equations above, keeping only linear terms 
in if and C- Of course, we must also linearize a e ff around 
its average (which is, in LD, just p): 



a e ff = p 



l eff 



at) 



(4) 



The result is a standard biased diffusion equation for ip 
within the lattice, along with novel equations for the ex- 
ceptional sites. In this approach, the diffusion coefficient 
is just 1/2 and the bias is v = 1 — 2p. Note that the only 
difference between /„ and I c lies in a^// m e 1- Q being 
zero or not. This set of linearized equations can be solved 
in Fourier space, with (p{k,Lo), ({lo), and ui). Defer- 
ring details to elsewhere, we present only the results along 
with a few remarks here. 

Assuming Gaussian noise with zero mean and corre- 
lation (£(fc,cj)£(fc ,u )} = AS k k >S u j , we solve the lin- 
earized equations and obtain, in particular, 



(ICMI 2 } = £ 



2A(1 - cosfc) 
|QHP(fc,u;)| 2 



where P(k, lu) = e lu — cos k + iv sin k and 

„ e luJ - 1 + a' ff (1 - p) (1 - e~ i<l ) 



At this level of approximation, A is usually taken as 
p (1 — p) [13] . Recalling that the complete power spectrum 
is precisely (|C(w)| 2 ), we can compare © with the simu- 
lation data. However, first we must account for a minor 
complication: the measurements taken only every I MCS. 
Summing over the unobserved t's is equivalent to summing 
over the harmonics w mi(tl = 2ir + 4) ; p — 0, ...,£ — 1. 
A straightforward computation leads to 



/(to) 



1 l ~ Y 

£2 



fi=0 k 



2A(l - cosfc) 



(7) 



(5) 



Note that the only difference between I Q and I c is the extra 
term proportional to ot' e j^ in Q(u>). 

In previous studies [8], a similar approach led to good 
fits, provided an effective "renormalized" diffusion coeffi- 
cient D is introduced and allowed as a fit parameter. Our 
result, as shown in eq. ©, appears more complicated than 
the previous approach. But, when it is evaluated numer- 
ically, it is very similar to the earlier prediction [8] for 
I (to). However, our UV cutoff is built in, so that there 
is no obvious way to insert a "renormalized" D into either 
eq. © or P(k,u). To go beyond this impasse, we plan to 
pursue perturbation theory seriously. Despite these unre- 
solved issues, let us present a remarkable prediction from 
eq. © here. Instead of the individual I's, consider the ra- 
tio Iq/I c - For reasons yet to be understood, there is quite 
good agreement between the data and theoretical results. 
In fig. [2l we illustrate this "zero parameter fit" for these 
ratios in a case with L — 1000. We should caution that 
the agreement in a larger system (L = 32000), though 
adequate, is not as impressive. Thus, the best we may 
conclude at this stage is that the ratio Iq/I c is somehow 
less sensitive to renormalization effects. 

Lastly, we turn to the crossover regime in the HD case, 
where a DW wanders relatively freely within the lattice. 
The appropriate comparison case in the ordinary TASEP 
is SP where I ~ to" 2 [8]. Our goal here is to explain the 
suppression of I c at small to. A physically intuitive picture 
is that the DW is localized by the feedback (through a e //) 
in the constrained TASEP [12]. Thus, its dynamics must 
be well described by a noisy harmonic oscillator and I c 
should be just a Lorentzian. For a quantitative picture, 
we start with eq. © for N p (t) = p(0,t), but with the 
simplifying approximations p(L) = 1-/3 and p(l) = a e //. 
For the SP in an ordinary TASEP, a e ff = a = /3, so 
that A t N p — £ (with appropriate boundary conditions), 
leading to I oc to~ 2 . With finite resources, a range of Ntot 
produces a e // [Np\ = f3, and so, we expand the second 
term to O (£)■ In the continuous t limit, eq. © becomes 
d t ( = -7C + C with 7 = (l-2j3)a' eff (N p ). The solution is 
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Fig. 2: (Color online) Simulation data and theory for the ra- 
tio of I and I c plotted as a function of m for LD (inset for 
crossover in HD). 



trivial, leading to (|£(u;)| 2 ) oc (w 2 + J 2 ) . As above, we 
compare the ratios I /I c , for which the theory predicts 



Io(m) 
I c (m) 



oc 1 + 



jT_ 

2-Km 



(8) 



In fig. [2l we plot this ratio from simulation data (for 
(L,a,/?,iV*,7V tot )=(1000,0.75,0.25,750,800)) as well as 
the expression above. The agreement is remarkably good, 
especially since no fit parameters (apart from the over- 
all proportionality constant) have been introduced. We 
have made such comparisons for other values of a, (3, etc. 
and all are similar, giving us confidence that this theory is 
quite adequate for predicting how the feedback mechanism 
suppresses I c (m) at small m. 

Summary and outlook. — We studied the power 
spectrum associated with the total occupation, N (t) , on 
a TASEP constrained by finite resources, using both simu- 
lations and analytic techniques. Considerable suppression 
at low frequencies is found, depending on the feedback 
through a e ff. Though much improved over a previous 
approach, the theory we studied also predicts only certain 
aspects of the observed spectra. Nevertheless, it is sur- 
prising that there is good agreement for the ratio of the 
power spectra (I /I c ), with no fit parameters! 

Obviously, this approach to understanding the simu- 
lation results leaves room for improvement. We believe 
that the non-linear terms neglected here will be the key 
to better predictions, especially at higher densities where 
the excluded volume constraint should be more relevant. 
Hopefully, a full investigation of their effects will also re- 
veal the secrets of the sensitive L-dependent of D and A 
found earlier [8, 14]. Beyond the study of a single TASEP 
coupled to a finite pool of particles, we look toward gen- 
eralizations of the model which may have further implica- 
tions for protein synthesis in a real cell. In particular, we 
should study TASEPs with large particles and inhomoge- 
neous hopping rates [3] . We are aware of preliminary stud- 
ies of the power spectra of these generalized TASEPs [15] 



and plan systematic investigations. Further, we should 
consider multiple TASEPs competing for the same pool 
of particles, just as many mRNA's in a cell compete for 
the same pool of ribosomes. A natural question - are there 
winners or losers? - can be crudely answered by studying 
their occupations, Ni (t), and currents, Ji (t). Obviously, 
power spectra of various combinations can reveal fluctu- 
ations and cross correlations between the different quan- 
tities. Perhaps they can be exploited as a sensitive di- 
agnostic tool for distinguishing different mechanisms that 
control translation. Beyond translation, we believe the 
study of power spectra will facilitate our understanding of 
non-equilibrium systems in other context, such as social 
networks, traffic flow, and finance. 
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